The Boids!¶

Flocking¶

The aggregate motion of a flock of birds, a herd of land animals, or a school of fish is a beautiful and familiar part of the natural world... The aggregate motion of the simulated flock is created by a distributed behavioral model much like that at work in a natural flock; the birds choose their own course. Each simulated bird is implemented as an independent actor that navigates according to its local perception of the dynamic environment, the laws of simulated physics that rule its motion, and a set of behaviors programmed into it... The aggregate motion of the simulated flock is the result of the dense interaction of the relatively simple behaviors of the individual simulated birds.

-- Craig W. Reynolds, "Flocks, Herds, and Schools: A Distributed Behavioral Model", Computer Graphics 21 4 1987, pp 25-34 See the original paper

Collision Avoidance: avoid collisions with nearby flockmates
Velocity Matching: attempt to match velocity with nearby flockmates
Flock Centering: attempt to stay close to nearby flockmates

Setting up the Boids¶

Our boids will each have an x velocity and a y velocity, and an x position and a y position.

We'll build this up in NumPy notation, and eventually, have an animated simulation of our flying boids.

In [1]:

import numpy as np

Let's start with simple flying in a straight line.

Our locations, for each of our N boids, will be an array, shape $2 \times N$, with the x positions in the first row, and y positions in the second row.

In [2]:

boid_count = 10

We'll want to be able to seed our Boids in a random position.

We'd better define the edges of our simulation area:

In [3]:

limits = np.array([2000,4000])

In [4]:

positions = np.random.rand(2, boid_count)*limits[:,np.newaxis]
positions

Out[4]:

array([[ 1676.65519964,  1582.67323371,   498.43795274,  1004.75928938,
         1253.99581583,   128.66738795,  1588.21375653,   106.9621584 ,
          129.25586159,   158.10522957],
       [ 1090.21820073,  3452.24003533,   657.02665987,   792.90879359,
         3482.15511352,  2197.79301965,  1569.78881754,  3120.12812311,
         3452.64033191,   581.66373878]])

In [5]:

positions.shape

Out[5]:

(2, 10)

We used broadcasting with np.newaxis to apply our upper limit to each boid. rand gives us a random number between 0 and 1. We multiply by our limits to get a number up to that limit.

In [6]:

limits[:,np.newaxis]

Out[6]:

array([[2000],
       [4000]])

In [7]:

limits[:, np.newaxis].shape

Out[7]:

(2, 1)

In [8]:

np.random.rand(2, boid_count).shape

Out[8]:

(2, 10)

So we multiply a $2\times1$ array by a $2 \times 10$ array -- and get a $2\times 10$ array.

Let's put that in a function:

In [9]:

def new_flock(count, lower_limits, upper_limits):
    width=upper_limits-lower_limits
    return (lower_limits[:,np.newaxis] + 
         np.random.rand(2, count)*width[:,np.newaxis])

But each bird will also need a starting velocity. Let's make these random too:

In [10]:

velocities = new_flock(boid_count, np.array([0, -20]), np.array([10, 20]))
velocities

Out[10]:

array([[  0.88195802,   1.98484018,   5.22024619,   7.03112532,
          9.76215136,   1.50556367,   9.55072034,   7.08702899,
          2.28081506,   7.12465118],
       [ 13.45373613, -12.9981505 ,   3.87008237,   4.01152776,
         -5.89533014,   6.37799663,  -9.27798715,   3.27817731,
        -17.4794183 ,   0.40560428]])

Flying in a Straight Line¶

Now we see the real amazingness of NumPy: if we want to move our whole flock according to

$\delta_x = \delta_t \cdot \frac{dv}{dt}$

We just do:

In [11]:

positions += velocities

Matplotlib Animations¶

So now we can animate our Boids using the matplotlib animation tools. All we have to do is import the relevant libraries:

In [12]:

from matplotlib import animation
from matplotlib import pyplot as plt
%matplotlib inline

Then, we make a static plot, showing our first frame:

In [13]:

# create a simple animation
positions=new_flock(100, np.array([100,900]), np.array([200,1100]))
velocities=new_flock(100, np.array([0,-20]), np.array([10,20]))

figure = plt.figure()
axes = plt.axes(xlim=(0, limits[0]), ylim=(0, limits[0]))
scatter=axes.scatter(positions[0,:],positions[1,:])
scatter

Out[13]:

<matplotlib.collections.PathCollection at 0x2afbc8e760f0>

Then, we define a function which updates the figure for each timestep

In [14]:

def update_boids(positions, velocities):
    positions += velocities
    
def animate(frame):
    update_boids(positions, velocities)
    scatter.set_offsets(positions.transpose())

Call FuncAnimation, and specify how many frames we want:

In [15]:

anim=animation.FuncAnimation(figure, animate,
                        frames=50, interval=50)

Save out the figure:

In [16]:

anim.save('boids_1.mp4')

And download the saved animation

You can even use an external library to view the results directly in the notebook. If you're on your own computer, you can download it from https://gist.github.com/gforsyth/188c32b6efe834337d8a (See the notes on installing libraries...)

Unfortunately, if you're on the teaching cluster, you won't be able to install it there.

In [17]:

from JSAnimation import IPython_display
# Inline animation tool; needs manual install via 
# If you don't have this, you need to save animations as MP4.
positions=new_flock(100, np.array([100,900]), np.array([200,1100]))
anim

Out[17]:

Fly towards the middle¶

Boids try to fly towards the middle:

In [18]:

positions=new_flock(4, np.array([100,900]), np.array([200,1100]))
velocities=new_flock(4, np.array([0,-20]), np.array([10,20]))

In [19]:

positions

Out[19]:

array([[  163.71294822,   194.34472245,   107.58100216,   198.50947897],
       [ 1012.44423944,   968.55940566,   964.95216618,   946.85987939]])

In [20]:

velocities

Out[20]:

array([[  1.63207235,   1.96849425,   0.8883695 ,   8.96225018],
       [  8.39194279,  -4.80451339, -18.52278931,  -9.03876946]])

In [21]:

middle=np.mean(positions, 1)
middle

Out[21]:

array([ 166.03703795,  973.20392267])

In [22]:

direction_to_middle = positions-middle[:, np.newaxis]

This is easier and faster than:

for bird in birds:
    for dimension in [0,1]:
        direction_to_middle[dimension][bird] = positions[dimension][bird]-middle[dimension]

In [23]:

move_to_middle_strength = 0.01
velocities = velocities - direction_to_middle * move_to_middle_strength

Let's update our function, and animate that:

In [24]:

def update_boids(positions, velocities):
    move_to_middle_strength = 0.1
    middle=np.mean(positions, 1)
    direction_to_middle = positions - middle[:, np.newaxis]
    velocities -= direction_to_middle * move_to_middle_strength
    positions += velocities

In [25]:

def animate(frame):
    update_boids(positions, velocities)
    scatter.set_offsets(positions.transpose())

In [26]:

anim=animation.FuncAnimation(figure, animate,
                        frames=50, interval=50)

In [27]:

from JSAnimation import IPython_display 
positions=new_flock(100, np.array([100,900]), np.array([200,1100]))
velocities=new_flock(100, np.array([0,-20]), np.array([10,20]))
anim

Out[27]:

Avoiding collisions¶

We'll want to add our other flocking rules to the behaviour of the Boids.

We'll need a matrix giving the distances between each bird. This should be $N \times N$.

In [28]:

positions=new_flock(4, np.array([100,900]), np.array([200,1100]))
velocities=new_flock(4, np.array([0,-20]), np.array([10,20]))

We might think that we need to do the X-distances and Y-distances separately:

In [29]:

xpos=positions[0,:]

In [30]:

xsep_matrix = xpos[:,np.newaxis] - xpos[np.newaxis,:]

In [31]:

xsep_matrix.shape

Out[31]:

(4, 4)

In [32]:

xsep_matrix

Out[32]:

array([[  0.        , -22.9607783 , -34.22480912, -66.73885357],
       [ 22.9607783 ,   0.        , -11.26403082, -43.77807527],
       [ 34.22480912,  11.26403082,   0.        , -32.51404445],
       [ 66.73885357,  43.77807527,  32.51404445,   0.        ]])

But in NumPy we can be cleverer than that, and make a 2 by N by N matrix of separations:

In [33]:

separations = positions[:,np.newaxis,:] - positions[:,:,np.newaxis]

In [34]:

separations.shape

Out[34]:

(2, 4, 4)

And then we can get the sum-of-squares $\delta_x^2 + \delta_y^2$ like this:

In [35]:

squared_displacements = separations * separations

In [36]:

square_distances = np.sum(squared_displacements, 0)

In [37]:

square_distances

Out[37]:

array([[    0.        ,  7216.85733881,  3048.95184722,  6984.18011626],
       [ 7216.85733881,     0.        ,  1605.96498346,  2908.15176693],
       [ 3048.95184722,  1605.96498346,     0.        ,  1105.72597211],
       [ 6984.18011626,  2908.15176693,  1105.72597211,     0.        ]])

Now we need to find birds that are too close:

In [38]:

alert_distance = 2000
close_birds = square_distances < alert_distance
close_birds

Out[38]:

array([[ True, False, False, False],
       [False,  True,  True, False],
       [False,  True,  True,  True],
       [False, False,  True,  True]], dtype=bool)

Find the direction distances only to those birds which are too close:

In [39]:

separations_if_close = np.copy(separations)
far_away = np.logical_not(close_birds)

In [40]:

separations_if_close[0,:,:][far_away] = 0
separations_if_close[1,:,:][far_away] = 0
separations_if_close

Out[40]:

array([[[  0.        ,   0.        ,   0.        ,   0.        ],
        [  0.        ,   0.        ,  11.26403082,   0.        ],
        [  0.        , -11.26403082,   0.        ,  32.51404445],
        [  0.        ,   0.        , -32.51404445,   0.        ]],

       [[  0.        ,   0.        ,   0.        ,   0.        ],
        [  0.        ,   0.        , -38.45889485,   0.        ],
        [  0.        ,  38.45889485,   0.        ,   6.96870761],
        [  0.        ,   0.        ,  -6.96870761,   0.        ]]])

And fly away from them:

In [41]:

velocities = velocities + np.sum(separations_if_close,2)

Now we can update our animation:

In [42]:

def update_boids(positions, velocities):
    move_to_middle_strength = 0.01
    middle = np.mean(positions, 1)
    direction_to_middle = positions - middle[:,np.newaxis]
    velocities -= direction_to_middle * move_to_middle_strength
    
    separations = positions[:,np.newaxis,:] - positions[:,:,np.newaxis]
    squared_displacements = separations * separations
    square_distances = np.sum(squared_displacements, 0)
    alert_distance = 100
    far_away = square_distances > alert_distance
    separations_if_close = np.copy(separations)
    separations_if_close[0,:,:][far_away] =0
    separations_if_close[1,:,:][far_away] =0
    velocities += np.sum(separations_if_close,1)
    
    positions += velocities

In [43]:

def animate(frame):
    update_boids(positions, velocities)
    scatter.set_offsets(positions.transpose())
    
anim = animation.FuncAnimation(figure, animate,
                        frames=50, interval=50)

from JSAnimation import IPython_display
positions = new_flock(100, np.array([100,900]), np.array([200,1100]))
velocities = new_flock(100, np.array([0,-20]), np.array([10,20]))
anim

Out[43]:

Match speed with nearby birds¶

This is pretty similar:

In [44]:

def update_boids(positions, velocities):
    move_to_middle_strength = 0.01
    middle = np.mean(positions, 1)
    direction_to_middle = positions-middle[:,np.newaxis]
    velocities -= direction_to_middle*move_to_middle_strength
    
    separations = positions[:,np.newaxis,:] - positions[:,:,np.newaxis]
    squared_displacements = separations*separations
    square_distances = np.sum(squared_displacements, 0)
    alert_distance = 100
    far_away=square_distances > alert_distance
    separations_if_close = np.copy(separations)
    separations_if_close[0,:,:][far_away] =0
    separations_if_close[1,:,:][far_away] =0
    velocities += np.sum(separations_if_close,1)
    
    velocity_differences = velocities[:,np.newaxis,:] - velocities[:,:,np.newaxis]
    formation_flying_distance = 10000
    formation_flying_strength = 0.125
    very_far=square_distances > formation_flying_distance
    velocity_differences_if_close = np.copy(velocity_differences)
    velocity_differences_if_close[0,:,:][very_far] =0
    velocity_differences_if_close[1,:,:][very_far] =0
    velocities -= np.mean(velocity_differences_if_close, 1) * formation_flying_strength
    
    positions += velocities

In [45]:

def animate(frame):
    update_boids(positions, velocities)
    scatter.set_offsets(positions.transpose())
    
anim=animation.FuncAnimation(figure, animate,
                        frames=200, interval=50)

from JSAnimation import IPython_display 

positions=new_flock(100, np.array([100,900]), np.array([200,1100]))
velocities=new_flock(100, np.array([0,-20]), np.array([10,20]))
anim

Out[45]:

Hopefully the power of NumPy should be pretty clear now. This would be enormously slower and, I think, harder to understand using traditional lists.